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Interactions between microorganisms and their complex flowing environments are essential in many 
biological systems. We develop a model for microswimmer dynamics in non-Newtonian Poiseuille 
flows. We predict that swimmers in shear-thickening (-thinning) fluids migrate upstream more (less) 
quickly than in Newtonian fluids and demonstrate that viscoelastic normal stress differences reorient 
swimmers causing them to migrate upstream at the centreline, in contrast to well-known boundary 
accumulation in quiescent Newtonian fluids. Based on these observations, we suggest a sorting 
mechanism to select microbes by swimming speed. 


Motile microorganisms ubiquitously inhabit confined 
and complex microenvironments. Geometrical con¬ 
straints are a key regulator of rheotaxis, the reorienta¬ 
tion of swimmers in response to external flows [2], and 
are essential in the design of microfluidic devices for drug 
delivery systems, hematology and cytometry [3, 4]. Ad¬ 
ditionally, the complexity of embedding fluids is crucial. 
One important aspect of the complexity arises from the 
dual fluidic and elastic (viscoelastic) behaviour of many 
biological fluids such as mucus and extracellular matrix 
gels [5-7] or blood at macroscopic length-scales [8-10]. 

Important correlations have been found between non- 
Newtonian behaviour of the fluid and pathological phe¬ 
nomena. Gastric mucus viscoelasticity effects swimming 
of H. pylori, an abundant pathogen in the stomach and 
leading cause of ulcers [11, 12]. It has been shown that 
viscoelasticity is a more crucial factor in controlling the 
maximum velocity of lyme disease pathogen B. burgdor¬ 
feri through skin than even chemical composition [13]. 
Viscoelastic properties of mucus have a remarkable im¬ 
pact on the swimming of spermatozoa and sperm-egg en¬ 
counter rates [14]. 

Despite the widespread implications of viscoelastic ef¬ 
fects on biological processes, research on motile microor¬ 
ganism dynamics in confined environments is largely lim¬ 
ited to Newtonian fluids [15-25]. Recently, a large num¬ 
ber of studies have considered locomotion in quiescent 
non-Newtonian fluids at the scale of microswimmers, in 
experiment, simulations and theory [26-33], but little is 
known about the dynamical behaviour of swimmers sub¬ 
ject to large-scale non-Newtonian flows. 

In this work, we construct a tractable theoretical 
framework for individual microorganisms swimming in 
confined, flowing microbiological environments of non- 
Newtonian fluids. We study the macroscopic effects 
of shear-dependent viscosity and viscoelasticity, both in 
separation and in conjunction, for a weakly viscoelas¬ 
tic fluid. Image systems are introduced, regularizing 
the hydrodynamic interaction of microswimmers with the 
walls, and swimmer trajectories are characterized. Shear- 
dependent viscosity is seen to greatly impact the up¬ 
stream motion of motile cells and our analysis shows that 
the presence of normal stress differences in viscoelastic 


fluids results in a remarkable upstream migration along 
the centreline. We provide quantitative measures of the 
upstream motion and propose a novel sorting mechanism 
for motile organisms in confined viscoelastic flows. 

A single microorganism is modeled as swimming in 
flowing, incompressible, non-Newtonian fluid within a 
channel of height 2H (Fig. 1). In addition to its swim¬ 
ming velocity in the direction Pg, the mo¬ 

tion of the swimming cell of radius a is affected by the 
background flow Vf, hydrodynamic interactions (HI) with 
the channel walls i;hi, and cross-streamline migration in¬ 
duced by viscoelastic normal stress difference gradients 
i;m- Thus, the evolution of a microswimmer’s position 
and direction are 

r^ = Vs+Vi + vhi + I’M (1) 

Pg = X Pg + X Pg, (2) 

where | V xvf and J7 hi denotes the angular velocity 

due to the HI with the walls. 

The translational invariance of Eqs. (1-2) along the y 
and z directions allows us to consider motion of swimmers 
in the ^ = 0 plane and orientation can be represented in 
cylindrical coordinates as Pg = — sin(0)ea, — co^{(j))ez^ 
where (f G [—7r,7r] is the angle in the x — z plane. 
Upstream swimming corresponds to 0 = 0 and down¬ 
stream to ±7r (Fig. 1). Gonsequently, the dynamics of 



FIG. 1. Schematic of a microswimmer at position Vs and mov¬ 
ing with speed Vs in the direction p^ subject to a viscoelastic 
flow within a microchannel of height 2H. The Poiseuille flow 
Vf is shown for shear-thinning (blue, dashed), Newtonian (red, 
solid) and shear-thickening (green, dotted) fluids. 
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the system can be represented by two coupled equations, 
X = x(x, (j)) and 0 0 ), and a third uncoupled equa¬ 

tion i = z{x^(j)). We nondimensionalise lengths by half 
the channel height, i7, and velocities by the swimming 
speed, v^. Therefore, changes in the swimming speed due 
to viscoelasticity, as studied in Refs. [26-33], are readily 
incorporated in this model. 

In a Newtonian fluid, this system shows the emergence 
of swinging and tumbling microswimmer trajectories in 
Poiseuille flow [17, 19]. Upstream-oriented swimmers are 
rotated by background vorticity so that they oscillate 
about the centreline (Fig. 2(a-b); green trajectory). For 
large oscillation amplitudes, however, the swimmer runs 
into the walls (Fig. 2(a); red trajectory). Hence, HI with 
the boundaries must be included [17]. Simply includ¬ 
ing the far-held force dipole of strength n and an image 
system consisting of a superposition of point-force singu¬ 
larities [34] in the HI produces non-physical singular how 
helds near the walls, unless a physical cut-oh length is 
provided. 

We construct a more physical representation by includ¬ 
ing a source doublet of strength a in the swimmer’s how 
and image helds, producing a more accurate near-held 
how and regularising the HI with the boundaries. This 
ensures that the swimmer is turned away from the bound¬ 
aries by the closest distance of approach 
which sets a natural cut-oh and gives an ehective size. 
This may be understood to be its hydrodynamic radius, 
ah = {2a [35], which we expect to be directly pro¬ 
portional to the swimmer size; ah ^ n, and thus hm ^ a. 
E.g. Volvox has a ^ 10^/rm^/s and Vs ^ lO^/rm/s [36], so 
that ah ^ 270/im compared to a ^ 200/im. By including 
the near-held correction, unphysical swimmer-wall con¬ 
tact is ruled out and the swimmer trajectory runs parallel 
to the wall with the ohset hm (Fig. 2(b); blue trajectory). 




(a) Non-regularised HI (b) Regularised HI 

FIG. 2. Typical trajectories for swimmer dynamics in a New¬ 
tonian Poiseuille how shown in x — 4> phase space, and in the 
X — z plane (insets). The swimmer and maximum how veloc¬ 
ities are I’s = 1 and I’max = 0.75, and the dipole moment is 
= 0. a) Non-regularised HI, cr = 0. b) Regularised HI with 
a — (3/10)^ so that hm = 3/10. The background colours 
indicate the velocity in the z direction. 


To consistently account for hnite size ehects, we also in¬ 
clude the Faxen corrections to the how induced transla¬ 
tional, Uf, and angular velocity, , of the swimmer [37] . 

Non-Newtonian ehects modify the background how 
and trajectories of microswimmers. Non-Newtonian hu- 
ids generally feature two properties diherent from a New¬ 
tonian counterpart — namely, shear dependent viscosity 
and normal stress diherences. Here shear-thinning and - 
thickening ehects are accounted for via a power-law huid 
model T] = '^ 70 ( 7 / 70 )’^”^: where 7 is the shear rate, 70 
is the viscosity at the shear-rate 70 , and n is the shear¬ 
thinning parameter. The background Poiseuille how of a 
power-law huid is 

t’f(r') = Vmax (1 - (|a;|/if)^) e^, (3) 

where 'r’max is the maximum how speed. This results in a 
stronger (weaker) how near the walls, in shear-thinning 
(-thickening) huids compared to a Newtonian huid with 
the same Vmax (Fig. 1). HI with the walls remain approxi¬ 
mately Newtonian for weakly non-Newtonian huids since 
the asymmetric correction for a dipolar swimmer [39, 40] 
decays rapidly as ^ r~^ [41, 42], which is small compared 
to the Newtonian contribution and amounts to a minor 
correction on the quadrupolar term. 

The upstream motion of a swimmer is enhanced in 
a shear-thickening huid compared to a shear-thinning 
counterpart without normal stresses (Fig. 3 and Supple¬ 
mental Material movie 1 [43]). This is associated with 
changes in vorticity in the vicinity of the walls. The 
stronger vorticity of the shear-thinning huid near the wall 
results in a more rapid reorientation towards the centre¬ 
line. Consequently, the swimmer has less time to move 
upstream. 

An initially upstream oriented swimmer (Fig. 3(a); 
blue trajectory) in a shear-tbinning huid moves a short 
distance upstream after the hrst oscillation about the 
centreline, whereas the swimmer in the shear-thickening 
huid progresses an order of magnitude further. Swim¬ 
mers initially orientated towards the walls (dashed 
green trajectories) are carried by the how, but in 
a shear-thickening huid they move further upstream 
near the walls. Similarly, swimmers initially orientated 
downstream (dotted red trajectories) experience an en¬ 
hanced downstream motion in a shear-t binning huid. 
This demonstrates that the dynamics in howing non- 
Newtonian environments can have a more signihcant ef¬ 
fect on motion than relatively small modihcations to the 
swimming speed in quiescent non-Newtonian huids [26- 
33]. 

If '^max = Vs, swimmers oriented directly upstream 
at the centreline do not progress, while those that os¬ 
cillate about the centreline experience less counterhow 
on average and therefore are able to migrate upstream 
(Fig. 3(b-c)). However, if the oscillations about the 
centreline are too large, the swimmer cannot move up¬ 
stream. Therefore, the ehective upstream motility is not 
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well described by any given trajectory but rather by a 
retention ratio [44], the ratio of the time-averaged 2 ;- 
component of swimmer velocity to the swimming speed, 
R = (t~^ Jq —z{t;xo^(t)o)dt^ where we average over 

all upstream-oriented trajectories xq G [—H hm^H — 
hm\ and G [—7r/2,7r/2]. The upstream retention ra¬ 
tio can be determined numerically (Fig. 3(d); points) 
and be approximated analytically. A conserved quan¬ 
tity of motion can be found by integrating x/0, giv¬ 
ing C = 1 + — VsCoscI). Hcncc, the dis¬ 

tance travelled along z per oscillation can be computed, 
D = idt, as well as the period, T = dt. Divid¬ 
ing these and averaging over the initial conditions gives 
R = — (:^) (Fig. 3(d); dashed lines). In the limit of 
'^max '^s we find the linear relation (Fig. 3(d); solid 
lines) 


2 2 + 9n -\- In? r’max 

TT 2-1- lOn + 12n^ 


( 4 ) 


Hence, the difference in upstream retention ratio for 



FIG. 3. Swimmer dynamics in Poiseuille flow of a shear¬ 
thinning (n = |) and -thickening fluid (n = 2) without nor¬ 
mal stresses, shown in the upper and lower halves of subfigures 
(a-c). Vs = r’max = 1, = 0 and a = (1/10)^. a) Trajectories 

in the x — z plane, with initial position Vg = (0,1/2) and orien¬ 
tations 0 = 0 (blue), 7r/2 (green), and tt (red), b) Trajectories 
in x — cf) phase space. The background colours indicate the ve¬ 
locity in the z direction, c) Upstream swimming velocity, —i, 
averaged over a large time, for all upstream-oriented initiali¬ 
sations in X — (j) space, d) Upstream retention ratio, defined 
by (c) averaged over these initial conditions, as a function of 
the flow speed. Points show full numerical solutions, dashed 
lines show theoretical predictions, and solid lines show the 
limit I’max ^ '^’s, Eq.(4). The inset focusses on this limit. 


shear-thinning and -thickening fluids grows with increas¬ 
ing flow speed. This determines the crossover between 
upstream or downstream motion of the majority of swim¬ 
mers where R = 0 (see inset of Fig. 3(d)). The slopes 
change at larger flow speeds, 'r’max > when the tum¬ 
bling trajectories start to outnumber the oscillating tra¬ 
jectories [17] and the full solution for R must be ap¬ 
plied (dashed lines). In this ^ Vs regime, the 

difference in upstream retention ratio for shear-thinning 
and -thickening fluids can be large (Fig. 3(d)). For 
'^max = lO'T’s, the shcar thickening (n = 2) i?-value differs 
by 33% from the shear-thinning (n = 1/2) value, which 
is substantial compared to the 5 — 10% change in swim¬ 
ming speed observed in quiescent non-Newtonian fluids 
[26-33]. 

The significant modification of upstream retention ra¬ 
tios in non-Newtonian fluids can have important conse¬ 
quences in microbiological flows. For instance, our re¬ 
sults suggest that a motile H. pylori^ swimming with an 
average velocity of 21[imjs [12] and subjected to gas¬ 
tric mucosal flow with a similar velocity and n = 0.5, 
would have a 50% reduction in upstream retention ra¬ 
tio than if it were swimming in a Newtonian fluid flow 
(n = 1). Since the velocity of the mucosal flow can vary 
broadly [45] and n can be as small as ^ 0.15 [12, 45, 46], 
this serves as a conservative example. 

In addition to shear-dependent viscosities, many mi¬ 
crobiological fluids are characterized by viscoelastic nor¬ 
mal stress differences. To describe these, with a power- 
law viscosity, we employ the second-order fluid model 
[47] with the stress tensor Sij = —pdij + 7^(7)D-j^ — 
-h (V^i + where t/i and 7/2 are the 

first and second normal stress coefficients and and 

( 2 ) 

D)d are the Rivlin-Eriksen tensors. The Deborah num- 
ber is De = ^ i xhe normal stress co- 

efficients characterize the fluid elasticity. These terms 
do not alter the undisturbed flow profile of Eq. (3) 
in the absence of swimmers. However, the disturbance 
flow around a finite-sized swimmer in combination with 
non-uniform shear across the channel results in a normal 
stress imbalance that causes a lateral migration across 
streamlines. Normal stress-induced migration of passive, 
inertialess particles in pressure-driven flow is well docu¬ 
mented [48-56]. To determine the migration velocity we 
use Chan and Leal’s solution for general quadratic flow 
[49] by expanding the background flow profile (Eq. 3) 
about the swimmer position, as reported previously [55] . 
In our system, the migration velocity is then 


Vm = -i^n 



( 5 ) 


where ”, /(n) = 5(1 + 

n)^“’^/36n^“’^ and ijjg = — 27 /^ 2 • The function 7 /^ 

encapsulates both the non-Newtonian effects of normal 
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stress differences and shear-dependent viscosity. A vis¬ 
coelastic torque flu is not included in Eq. (2) [55] be¬ 
cause this term is not significant compared to the vor- 
ticity when De ^ 1 and, by symmetry of the swimmer, 
does not lead to preferred orientations. 

In both the shear-thinning and -thickening cases with 
normal stresses, the swimmer is driven to the centreline, 
and the coupling between motility and streamline migra¬ 
tion rotates the swimmer to move upstream along the 
centreline (Fig. 4(a)). Unlike in a Newtonian fluid, the 
oscillations about the centreline are now damped in am¬ 
plitude as the phase space origin (x = 0 = 0) is a stable, 
attractive spiral (Fig. 4(b)). The attraction is stronger 
for shear-thinning than shear-thickening fluids. 

We analyse this effect by linearising the equations of 
motion (1-2) about the origin so that (0, x)^ = x)^ 

where 


M = 


■ 4 


V{+^ 


_L ^ I ^ — 7/; 

^3141 2 2 


(6) 


In AT, and dipolar HI terms are responsible for the 
spiral. Away from the walls, viscoelasticity dominates 
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FIG. 4. Swimmer dynamics in Poiseuille flow of a shear¬ 
thinning viscoelastic fluid with De = 0.1, n — 0.8, and 
r’max = Ir’s- The swimmer parameters are a = 0.1, Vs = 1, 
= 0 and a = 10~^. a) Oscillating trajectory in the 
X — z plane with initial position rg = (0,0.9) and orienta¬ 
tions 0 = 0. b) Corresponding trajectory in x — 0 phase 
space. Colours indicate time progressing, from t = 0 (blue) 
to t = 1000 (red). The swimmer is focussed towards the 
centerline and are reoriented to move upstream, c-d) Two 
ensembles of swimmers, with Vs = 1.1 (blue) and Vs = 0.9 
(red), are released from a random x-position and orientation 
in the channel at z = 0. In a Newtonian fluid (c), the swim¬ 
mers are dispersed but in a viscoelastic fluid (d) they remain 
clustered and are sorted according to swimming speed over 
time. 


over HI effects and the eigenvalues of M without HI are 
found to be X± = ± 4'^max'^s)- Hence, the 

origin is a stable fixed point if 7 /;^ > 4i;max'^s with two real 
negative (attractive) eigenvalues. Otherwise, the origin 
is a stable spiral with complex eigenvalues and negative 
real parts, meaning that swimmers perform damped os¬ 
cillations about the centreline as verified in Fig. 4(a,b). 
Because the function f{n) decreases monotonically with 
n, ipn is larger for shear-thinning fluids and therefore the 
attraction towards the centreline is greater. 

Though more pronounced in shear-thinning than 
shear-thickening flows, swimmers in flowing viscoelas¬ 
tic fluids tend to move upstream along the centreline 
after some time, regardless of initial position or up¬ 
stream orientation. This allows for a sorting mecha¬ 
nism to select swimmers with a given swimming speed 
larger than the tunable Poiseuille flow, as demonstrated 
in Fig. 4(c,d), where distributions of swimmers with dif¬ 
ferent self-propulsion velocities are initially introduced 
at random positions and orientations in the channel in 
Newtonian and shear-thinning viscoelastic fluids. Unlike 
the Newtonian fluid, swimmers with larger motility are 
separated by moving upstream in the viscoelastic fluid 
(see Supplemental Material movies 2,3 in Ref. [43]). It is 
worth noting that we expect this sorting mechanism to 
be robust against translational and orientational noise 
since small amounts of noise will keep the oscillation 
size nonzero, enhancing the upstream retention ratio and 
hence the sorting. 

To summarize, unlike the prevalent boundary accumu¬ 
lation in quiescent Newtonian fluids, swimmers’ trajec¬ 
tories show oscillatory motion about the centreline. Av¬ 
erage migration against Poiseuille flows is enhanced (re¬ 
duced) in shear-thickening (-thinning) fluids compared 
to simple Newtonian fluids. It is not necessary that the 
non-Newtonian nature of these fluids be appreciable on 
the microscale since altered trajectories arise from differ¬ 
ences in vorticity at macroscopic scales. This constitutes 
a substantial change to the effective upstream motility, 
that is comparable to or greater than observed changes 
in motility due to microscopic effects on swimming in 
quiescent non-Newtonian fluids [26-33]. The oscillations 
are damped towards the centreline in the presence of vis¬ 
coelastic normal stress differences resulting in direct up¬ 
stream migration. This offers a sorting mechanism to 
differentiate motile microorganisms according to speed. 
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